%% it computes the invariant distribution of b

b_vec_l = (log(0.01):db:log(1+db))+log(b_bar);
b_vec_d = exp(b_vec_l);
N=length(b_vec_d);

h_00     = zeros(N,1);
inde_hig = find(b_vec_d>bhig,1);
inde_low = find(b_vec_d>blow,1);
inde_mid = find(b_vec_d>bmid,1);
h_00(inde_hig)=q_hig*incid;
h_00(inde_low)=q_low*incid;
h_00(inde_mid)=q_mid*incid;
if incid<1
    inde_hig = find(b_vec_d>b_up2,1);
    inde_low = find(b_vec_d>b_dn2,1);
    inde_mid = find(b_vec_d>b_md2,1);
    h_00(inde_hig)  =h_00(inde_hig)+q_hig*(1-incid);
    h_00(inde_low)  =h_00(inde_low)+q_low*(1-incid);
    h_00(inde_mid)  =h_00(inde_mid)+q_mid*(1-incid);
end

%%
Tb_f = zeros(N,N);
for b=1:N-1
    Tb_f(b,b)   = -1/db;
    Tb_f(b,b+1) =  1/db;
end
Tb_f=sparse(Tb_f);

Tb_b = zeros(N,N);
for b=1:N-1
    Tb_b(b+1,b+1)=   1/db;
    Tb_b(b+1,b)  =  -1/db;
end
Tb_b = sparse(Tb_b);
Tbb  = zeros(N,N);
Tbb(1,1)      =  -1/db^2;
Tbb(1,2)      =   1/db^2;
Tbb(N,N)      =  -1/db^2;
Tbb(N,N-1)    =   1/db^2;
for b=2:N-1
    Tbb(b,b-1)=  1/db^2;
    Tbb(b,b)  = -2/db^2;
    Tbb(b,b+1)=  1/db^2;
end
Tbb=sparse(Tbb);

drift_I   =  ppval(pp_drift,b_vec_d);

if inde_dilution ==1
    ell        = ppval(pp_ell,b_vec_d) ;ell(end)=ell(end-1);
else
    ell        = zeros(1,length(b_vec_d));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
inde_b     = (ell./b_vec_d-(rho+drift_I-sigma^2/2))';inde_b(1)=inde_b(2);inde_b(end)=inde_b(end-1);
Tb         = Tb_f;
Tb_db      = sparse(Tb.*(inde_b*ones(1,N)));
Tbb_dbb    = sparse(Tbb*(sigma^2/2));

inde_bb    = find(b_vec_d>b_bar);
inde_bi    = find(b_vec_d<=b_bar);
inde_nego  = find(b_vec_d>=b_bar*a,1);
A_exit     = speye(N)*delta;
AA         = Tb_db + Tbb_dbb;
E = zeros(1,N);
E(1,inde_nego)=1;
EE = (repmat(sparse(E),N,1));%matrix of zeros, except for the location of the renegotiation
X  = sum(AA(:,inde_bb),2);%vector of mass of exiters from each greed point j
X(inde_bb)=0;
XX  = repmat(X,1,N);%repeat the vector in multiple equal columns
HH  = (repmat(sparse(h_00'),N,1));

%obtain the transition matrix with renegotiation
AA2 = -A_exit + AA + phi*(EE.*XX);%here add mass only when XX(j,i) is positive, which is when there is exit from point (j) AND EE(j,i)=1 that is when i is a renegotiation point given j
AA2(inde_bb,inde_bb)=-speye(length(inde_bb));

%transpose the transition matrix so that you go from HJB to KFE; the minus sign is to take it to left of the equation
A_T = -AA2';
%% find the distribution
crit_h=1;sum_h=1;
h_til     = zeros(N,1);
while crit_h>1/1000000000000
    h_til(inde_bi)  = A_T(inde_bi,inde_bi)\(h_00(inde_bi)*(1+sum_h*delta));
    crit_h = max(abs(sum(h_til)-sum_h));
    sum_h=sum(h_til);
end

bnkptcy_rate  = 1/sum_h;
exit_rate     = bnkptcy_rate + delta;

f_dist = h_til(inde_bi);
b_dist = b_vec_d(inde_bi);

